# aem.two_d.transient.head_line_sink -- 
#
# Copyright (c) 2000, V.Kelson and M.Bakker   (steady head line sink)
#           (c) 2001, Willem Jan Zaadnoordijk (transient head line sink)

from types import *

import geometry
XY = geometry.XY
SearchableList = geometry.two_d.collections.SearchableList

import base

import aem.container
RowList = aem.container.RowList
ColumnList = aem.container.ColumnList

import aem.matrix



class HeadLineSink(geometry.two_d.Line,base.Element):
	'''Class for transient constant head linesinks

	Python call: HeadLineSink(pt1,pt2,tHlist=[geometry.XY(0,0)],parent=None)

	Attributes:

	  - pt1: startpoint coordinates
	  - pt2: end  point coordinates
	  - width: width of line (from Line)
	  - z1: complex coordinate of pt1 (from Line)
	  - z2: complex coordinate of pt2 (from Line)
	  - zc: complex coordinate of center of line (from Line)
	  - zl: complex coordinate 0.5(z2-z1) (from Line)
	  - halfLength: half the length of the line (from Line)
	  - parameters: ColumnList of strength parameters, set to empty list (inherited from Element)
	  - solver: solver (inherited from Element)
	  - tHlist:  list of pairs (time, head) with time in ascending order
	  - parent: parent of the element (inherited from Element)
	  - elementList: list of transient elements ("UnitStepLineSink"s)
	  - initialSteadyState: list of initial steady element(s)
	    this element is not included in "parameters"and "elementList"
	    and does not appear in the influence functions

	Overloaded functions from Element:
		__init__, getMatrixCoefficients, getMatrixContribution,
		check, potentialInfluence, dischargeInfluence
	'''

	def __init__(self,pt1,pt2,tHlist=[geometry.XY(0,0)],parent=None,overspec=1):
		'''Constructor'''
		geometry.two_d.Line.__init__(self,pt1,pt2)
		base.Element.__init__(self,parent)
		self.elementList = SearchableList()
		self.__numberOfParameters = 0
		self.parameters= ColumnList([])
		self.initialH  = None
		self.tlastsolve= -1
		self.specifiedH= []
		self.specifiedt= []
		self.givenElement = 0
		# Decode options for control-points  #function inherited from Line
		self.controlPoints = self.computeControlPoints(overspec)
		for partim in tHlist :
			if partim.x <= 0 :
				if self.tlastsolve < 0 :
					# no initial steady state has been specified yet
					self.initialH = partim.y
					self.initialSteadyState = \
						aem.two_d.single.HeadLineSink(pt1,pt2,partim.y,
						                      parent.initialSteadyState)
					# initial steady state is starting head for transient calculation
					self.tlastsolve = 0
					self.specifiedH.append(partim.y)
					self.specifiedt.append(partim.x)
			else :
				self.specifiedH.append(partim.y)
				self.specifiedt.append(partim.x)
				if self.tlastsolve < 0 : selftlastsolve = partim.x

	def __repr__(self):
		return 'Trans.Head LineSink{'    \
		   + str((self.pt1,self.pt2,     \
		         self.initialH,self.specifiedt,self.specifiedH)) \
		   + '}'

	def addElement(self,el):
		self.elementList.append(el)

	def getMatrixCoefficients(self,fa,ctlPt,t):
		# return value of the FA influence function at point CTLpT and time T
		# for the last UnitStepLineSink of this element

		if t > self.tlastsolve :
			for i in range(len(self.specifiedt)) :
				if t > self.specifiedt[i] :
					# determine start time connected to this solve time and
					# add UnitStepLineSink if necessary
					# - same code as in getMatrixContribution, because
					#   it is unknown which function gets called first
					nelem=len(self.elementList)
					if nelem==0:
						tlast=self.tlastsolve-1
					else:
						elemlast = self.elementList[len(self.elementList)-1]
						tlast=elemlast.t
					if tlast<self.tlastsolve:
						tstart = max(self.tlastsolve,self.specifiedt[0])
						aem.two_d.transient.UnitStepLineSink(         \
							self.halfLength,tstart,self)
						self.parameters.append(0)
						self.__numberOfParameters = self.__numberOfParameters+1
					# start time has been determined (and UnitStepLineSink added)
					# - same code as in getMatrixCoefficients, because
					#   it is unknown which function gets called first
		cz=complex(ctlPt.x,ctlPt.y)
		czl = (2*cz-(self.z2+self.z1))*self.halfLength / (self.z2-self.z1)
		ctlPtlc = XY( czl.real , czl.imag )
		if len(self.elementList)>0:
			return fa(self.elementList[(len(self.elementList)-1)],ctlPtlc,t)
		else:
			return 0

	def getMatrixContribution(self,elementList,t):
		# Returns a tuple (mc,mr) where mc is a list of MatrixColumn objects
		# and mr is a list of MatrixRow objects
		# specifying the condition at time t of a new
		# UnitStepLineSink starting at the last time of solving
		rows = []
		columns = []
		specifiedHead= None
		active = 0
		if t > self.tlastsolve :
			for i in range(len(self.specifiedt)) :
				if t > self.specifiedt[i] :
					specifiedHead = self.specifiedH[i]
					active = 1
		if active==1 :
			# determine start time connected to this solve time and
			# add UnitStepLineSink if necessary
			# - same code as in getMatrixContribution, because
			#   it is unknown which function gets called first
			nelem=len(self.elementList)
			if nelem==0:
				tlast=self.tlastsolve-1
			else:
				elemlast = self.elementList[len(self.elementList)-1]
				tlast=elemlast.t
			if tlast<self.tlastsolve:
				aem.two_d.transient.UnitStepLineSink(         \
					self.halfLength,self.tlastsolve,self)
				self.parameters.append(0)
				self.__numberOfParameters = self.__numberOfParameters+1
			# start time has been determined (and UnitStepLineSink added)
			# - same code as in getMatrixContribution, because
			#   it is unknown which function gets called first
			columns.append(aem.matrix.MatrixColumn(self,len(self.elementList)-1))
			for cp in self.controlPoints:
				self.alpha=self.parent.aquiferDiffusivity(cp,t)
				mr = aem.matrix.MatrixRow()
				for e in elementList:
					mr.A = mr.A + e.getMatrixCoefficients(lambda z,p,t:z.potentialInfluence(p,t),cp,t)
					mr.b = self.headToPotential(cp,specifiedHead) - \
						            self.potential(cp,t)
			rows.append(mr)
		self.tlastsolve = t
		return ( rows,columns )

	def check(self,t):
		# Checks the boundary condition(s)
		rv = []
		active = 0
		for i in range(len(self.specifiedt)) :
			if t > self.specifiedt[i] :
				specifiedHead= self.specifiedH[i]
				active = 1
		if active==1 :
			for cp in self.controlPoints:
				rv.append(self.head(cp,t) - specifiedHead)
		return rv

	def potentialInfluence(self,xy,t):
		assert isinstance(xy,XY),"XY point is required"
		pot=[]
		self.alpha=self.parent.aquiferDiffusivity(xy,t)
		cz=complex(xy.x,xy.y)
		czl = (2*cz-(self.z2+self.z1))*self.halfLength / (self.z2-self.z1)
		xyw = XY( czl.real , czl.imag )
		for el in self.elementList:
			pot = pot + el.potentialInfluence(xyw,t)
		return RowList( pot )

	def dischargeInfluence(self,xy,t):
		assert isinstance(xy,XY),"XY point is required"
		dis=[]
		self.alpha=self.parent.aquiferDiffusivity(xy,t)
		cz=complex(xy.x,xy.y)
		czl = (2*cz-(self.z2+self.z1))*self.halfLength / (self.z2-self.z1)
		xyw = XY( czl.real , czl.imag )
		for el in self.elementList:
			disell= el.dischargeInfluence(xyw,t)
			for d in disell :
				rqxi = d.x
				rqeta= d.y
				rqx = ( (self.z2-self.z1).real*rqxi   \
				       -(self.z2-self.z1).imag*rqeta ) / abs(self.z2-self.z1)
				rqy = ( (self.z2-self.z1).imag*rqxi   \
				       +(self.z2-self.z1).real*rqeta ) / abs(self.z2-self.z1)
				d = XY(rqx,rqy)
			dis = dis+disell
		return RowList( dis )
